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Abstract. The convergence and optimality of adaptive mixed finite element methods 
for the Poisson equation are established in this paper. The main difficulty for mixed finite 
element methods is the lack of minimization principle and thus the failure of orthogonal- 
ity. A quasi-orthogonality property is proved using the fact that the error is orthogonal 
to the divergence free subspace, while the part of the error that is not divergence free can 
be bounded by the data oscillation using a discrete stability result. This discrete stability 
result is also used to get a localized discrete upper bound which is crucial for the proof 
of the optimality of the adaptive approximation. 



Adaptive methods are now widely used in scientific computation to achieve better 
accuracy with minimum degrees of freedom. While these methods have been shown to be 
very successful, the theory ensuring the convergence of the algorithm and the advantages 
over non-adaptive methods is still under development. Recently, several results have been 
obtained for standard finite element methods for elliptic partial differential equations 
[8, 36, 48,50, 12, 59,51,27,29]. 
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In this paper, we shall establish the convergence and optimality of adaptive mixed 
finite element methods (AMFEMs) of the model problem 

- Am = / in f2, and u = on <9f2, (1.1) 

posed on a polygonal and simply connected domain f2 C IR 2 . In many applications ([24]) 
the variable cr = — Vu is of interest and it is therefore convenient to use mixed finite 
element methods, such as the Raviart-Thomas mixed method [53] and Brezzi-Douglas- 
Marini mixed method [23]. We shall construct adaptive mixed finite element methods 
based on the local refinement of triangulations and prove they will produce a sequence 
of approximation of cr in an optimal way. 

Our main result is the following optimal convergence of our algorithms AMFEM and 
its variant. Let cr N be the approximation of cr based on the triangulation Tn obtained in 
AMFEM. If cr G A 8 and / G A S Q , then 

Ik - <t n \\ < C{\\<t\\ a . + ||/|U g )(#7]v - #T )- S , (1.2) 

where (A s , \\ ■ and (A S Q , \\ ■ \\a») are approximation spaces as in [12]. The index 
s is used to characterize the best possible approximation rate of cr, which depends on 
the regularity of the solution and data, and the order of elements. For example when 
/ G L 2 (Vt) and cr G W l,1 (VL), we can achieve the optimal convergence rate s = 1/2 for 
the lowest order Raviart-Thomas finite element space. We refer to [13] for the charac- 
terization of A s in terms of Besov spaces and to [9, 10, 35, 34] for the regularity results 
in Besov norms. We comment that to apply our adaptive algorithm, we do not need to 
know s explicitly. Our algorithm will produce the best possible approximation rate for 
the unknown cr. 

For the analysis of the convergence of adaptive procedure, we follow the new approach 
by Cascon, Kreuzer, Nochetto and Siebert [27], and for the optimality we mainly use the 
simplified case in Stevenson's work [59]. A distinguish feature of the new approach 
for the convergence proof is the relaxation of the interior node requirement for the re- 
finement. We do not claim any originality on the proof of convergence and optimality. 
Instead the main contribution of this paper is to establish two important ingredients used 
in the proof, namely quasi-orthogonality and discrete upper bound. 

One main ingredient in the convergence analysis of standard AFEM is that the error 
is orthogonal to the finite element spaces in energy-related inner product since the stan- 
dard finite element approximation can be characterized as a minimizer of Dirichlet-type 
energy. For mixed finite element methods, however, the approximation is a saddle point 
of the corresponding energy and thus there is no orthogonality available. We shall prove 
a quasi-orthogonality result. A similar result for the lowest order Raviart-Thomas finite 
element space has recently been proved by Carstensen and Hoppe [26], where a special 
relation between mixed finite element method and non-conforming method is used. In 
this paper, we shall propose a new and more straight-forward approach which works for 
any order elements and both Raviart-Thomas and Brezzi-Douglas-Marini methods. The 
main observation is that the error is orthogonal to the divergence free subspace, while 
the part of the error containing divergence can be bounded by the data oscillation using 
a discrete stability result. 

Another ingredient to establish the optimality of the adaptive algorithm is the localized 
discrete upper bound for a posteriori error estimator. Using the discrete stability result, 
we are able to obtain such discrete upper bound and use it to prove the optimality of the 
convergent algorithm. The optimality of mixed adaptive finite element methods seems to 
be new. 
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The rest of this paper is organized as follows. In Section 2, we shall introduce mixed 
finite element methods and give a short review of mesh adaptivity through local refine- 
ment. We shall include many preliminary results in this section for later usage. In Section 
3, we shall prove the discrete stability result and use it to prove the quasi-orthogonality 
result. In Section 4, we shall present a posteriori error estimator and prove the discrete 
upper bound. In Section 5, we shall present our algorithms and prove their convergence 
and optimality. 

Throughout this paper, we shall use standard notation for Sobolev spaces and use 
boldface letter for the spaces of vectors. The letter C, without subscript, denotes generic 
constants that may not be the same at different occurrences and Cj, with subscript, de- 
notes specific important constants. 

2. Preliminaries 

In this section we shall introduce mixed finite element methods for the Poisson equa- 
tion and discuss the general procedure of adaptive methods through local refinement. We 
shall also include a result on the approximation of the data. 

2.1. Mixed finite element methods. The standard finite element method involves writ- 
ing (1.1) as a primal variational formulation: for a given / G L 2 (Vl), find u G Hq(Q) 
such that 

/ V« • Vv = [ fv Vv G H^(n), (2.1) 
Jn Jn 

and then finding an approximation by solving (2.1) in finite-dimensional subspaces of 
Hq(CI). In many applications ([24]) the variable cr = — Vw is of interest, and it is 
therefore convenient to use mixed finite element methods. Let us first write (1.1) as a 
first order system: 

cr + Vu = 0, div cr = f mil, and u = on Oil. (2.2) 

Let 

£ = JET(div;ft) := {r G L 2 (tt) : divr G L 2 (tt)}, and U = L 2 (VL). 
We shall use || • || to denote L 2 -norm and || • ||_f/(div) f° r the Jf(div) norm: 

||T||^ (div ) = (||T|| 2 +||divT|| 2 ) 1 / 2 , VrGE. 

The mixed (or dual) variational formulation of (2.2) is, given an / G L 2 {VL), find 
{cr,u) G £ x U such that 

(cr, t) - (div r, u) = Vr G S, (2.3) 

(div<r,u) = (f,v) VveU, (2.4) 

where (-, ■) is the inner product for L 2 {Vt) or L 2 (Q). Note that the Dirichlet boundary 
condition is imposed as a natural boundary condition in the dual formulation (2.3) using 
integration by parts. The existence and uniqueness of the solution (cr,u) to (2.3)-(2.4) 
follows from the so-called inf-sup condition which can be easily established for this 
model problem [24]. 

Given a shape regular and conforming (in the sense of [30]) triangulation Th of fi, the 
mixed finite element method is to solve (2.3)-(2.4) in a pair of finite-dimensional spaces 
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Ti H C £ and U H C U. That is, given an / G L 2 (tt), to find (cr H ,u H ) G T, H x [/# such 
that 

(<t h , t h ) - (div t h , u H ) = Vt h 6 S h (2.5) 

(div a H , v H ) = Uh,v h ) V% 6 (2.6) 

Hereafter /# denotes the L 2 (Q) projection of / onto U H . Namely, fn^Un such that 
(/hjVh) = (/,Vh), Vuh G Uh- The well-posedness of the discrete problem (2. 5)-(2. 6), 
unlike the standard finite element method for the primary variational formulation, is non- 
trivial. One sufficient condition to construct stable finite element spaces is to ensure the 
inf-sup condition still holds for the discrete problem. Since 1970's many stable finite 
element spaces have been introduced for this case, such as those of Raviart-Thomas 
spaces [53] and Brezzi-Douglas-Marini spaces [23]. Recently it has been shown that 
such stable finite element spaces can be constructed in an elegant way using differential 
complex theory [16, 41, 2, 5]. 

The Raviart-Thomas spaces [53] are defined for an integer p > by 

RT H = YP H x U*, where 
S^(Th) := {r G H(dw;Q) : r\ T G P p (T) + xV p (T), VT G T H }, 
and U P H {T H ) := {v G L 2 (Q) : v\ T G V P (T), VT G T H }, 

and where V P (T) denotes the space of polynomials on T of degree at most p. 
The Brezzi-Douglas-Marini spaces [23] are defined for an integer p > 1 by 

BDM H = 'S p H xUf I , where 

S^(Th) := {r G #(div;Q) : r\ T G P p (T), VT G Th}, 
and U*{T H ) ■= {v G L 2 (Q) : v\ T G P p _i(T), VT G 7^}. 

Since most results hold for both Raviart-Thomas and Brezzi-Douglas-Marini spaces 
and p is fixed in most places, we shall use generic notation (E#, Uh) to denote the pair 
in RT H or BDM H - The discrete problem posed on (E#, C/^) will satisfy the discrete 
inf-sup condition [24] from which the existence and uniqueness of the finite element 
approximation (cth,uh) follows. 

We shall use £ and Ch to denote the differential operators corresponding to (2.3)-(2.4) 
and (2.5)-(2.6), respectively. Those equations can be formally written as 

£(<t,u) = f and C h ((t h ,u h ) = f H . 

We shall use the notation (cr,u) = £ -1 / and (cth,uh) = C>h Ih to emphasis the 
dependence of /. With an abuse of notation, we also use cr = C~ l f and cr H — ^h 
when cr and a H are of interest. 

2.2. Adaptive methods through local refinement. Let cr = Cr x j and <r H = fn- 

We are mostly interested in the control of the error \\<r — <th\\ which is usually more 
important than control the error of scalar variable u in mixed finite element methods. 
Although the natural norm for the error is \\cr — er#||#( div ), we comment that, by (2.4) 
and(2.6), || diver — div cr H || = ||/ — can be approximated efficiently without solving 
equations and also may dominate the error ||er — cr#||#( div ); see Remark 3.4 in [45]. 

The rate of the error ||<r — <th\\ for cth G H p h (Th) depends on the regularity of the 
function being approximated and the regularity of the mesh. If cr G i? p+1 (fi) and Th is 
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quasi-uniform with mesh size H = max^-j^ diam(T), then the following convergence 
result of optimal order is well known [24] 

Ik - tr H \\ < CH p+1 \\a\\ p+1 . (2.7) 

The regularity result cr e H P+1 (Q), however, may not be true in many applications, 
especially for concave domains Q. Thus we cannot expect the convergence result (2.7) 
on quasi-uniform grids in general. 

To improve the convergence rate, element sizes are adapted according to the behavior 
of the solution. In this case, the element size in areas of the domain where the solution 
is smooth can stay bounded well away from zero, and thus the global element size is 
not a good measure of the approximation rate. For this reason, when the optimality of 
the convergence rate is concerned, #7~, the number of elements, is used to measure the 
approximation rate in the setting of adaptive methods that involve local refinement. 

We now briefly review the standard adaptive procedure. Given an initial triangula- 
tion %, we shall generate a sequence of nested conforming triangulations Tk using the 
following loop 

SOLVE -> ESTIMATE -> MARK ->■ REFINE. (2.8) 

More precisely, to get Tk+i from Tk we first solve (2.5)-(2.6) to get cr fc on Tk- The error 
is estimated using er fc and data. And the error estimator is used to mark a set of of 
triangles or edges that are to be refined. Triangles are then refined in such a way that the 
triangulation is still shape regular and conforming in the sense of [30]. 

We shall not discuss the step SOLVE which deserves a separate investigation. We 
assume that the solutions of the finite-dimensional problems can be generated to any ac- 
curacy to accomplish this in optimal space and time complexity. Multigrid-like methods 
for mixed finite element methods on quasi-uniform grids can be found in [17, 18, 20, 21, 
40,52,56]. 

The a posteriori error estimators are essential part of the ESTIMATE step. Given a 
shape regular triangulation Th, let Eh denote the edges of Th ■ In this paper, we shall use 
edge-wise error estimator t)e for each edge E E Eh- See Section 4 for details. 

The local error estimator t)e is employed to mark for refinement the elements whose 
error estimator is large. The way we mark these triangles influences the efficiency of the 
adaptive algorithm. In the MARK step we shall always use the marking strategy firstly 
proposed by Dorfler [36] in order to prove the convergence and the optimality of the local 
refinement strategy. 

In the REFINE step we need to carefully choose the rule for dividing the marked tri- 
angles such that the mesh obtained by this dividing rule is still conforming and shape 
regular. Such refinement rules include red and green refinement [11], longest edge re- 
finement [55, 54], and newest vertex bisection [58, 46, 47]. Note that not only marked 
triangles get refined but also additional triangles are refined to recovery the conformity 
of triangulations. We would like to control the number of elements added to ensure the 
overall optimality of the refinement procedure. To this end, we shall use newest vertex 
bisection in this article. We refer to [46, 61, 12, 28] for details of newest vertex bisection 
and only list two important properties below. 

Let % be a conforming triangulation refined from a shape regular triangulation % 
using the new vertex bisection and let Ai be the collection of all marked triangles going 
from % to Tk. Then 

(1) {Tk} is shape regular and the shape regularity only depends on %; 

(2) + C#M. 
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Recently Stevenson [60] showed that such results can be extended to bisection algo- 
rithms of n-simplices. The optimality of the adaptive finite element method in this paper, 
thus, could be extended to general space dimensions. 

2.3. Approximation of the data. We shall introduce the concept of data oscillation 
which is firstly introduced in [48], and use it here for the approximation of data. Such 
quantity measures intrinsic information missing in the averaging process associated with 
finite elements, which fails to detect fine structures of /. 

For a set A, Ha denotes the diameter of A. To simplify the notation, we may drop 
the subscript if it is clear from the context. For a triangulation Th of Vt and a function 
/ G L 2 (n), we define a triangulation dependent norm 

ll^/llo,rH:=(E^II/llo,T) 1/2 - 

T£T H 

Definition 2.1. Given a shape regular triangulation Th of Q and an f G L 2 (Q), we 
define the data oscillation 

osc(f,T H ):= \\H(f - f H )\\ 0)TH . 

Let Vn denote the set of triangulations constructed from an initial triangulation To by 
the newest vertex bisection method with at most N triangles. We define 




where iV is a fixed integer representing the number of triangles in %. We will recall a 
result of Binev, Dahmen and DeVore [12] which shows that the approximation of data 
can be done in an optimal way. The proof can be found at [12]; See also [14]. 

Theorem 2.2 (Binev, Dahmen and DeVore). Given a tolerance e, an f G L 2 (Q) and a 
shape regular triangulation To, there exists an algorithm 

T H = APPROX(/,T ,£) 

such that 

osc(/,Th)<£, and #T H -#%<C\\ffj, s e- 1/s . 

3. QUASI-ORTHOGONALITY 

Unlike the primal formulation of Poisson equation, cr H is not the L 2 -orthogonal pro- 
jection of cr from £ to Indeed the solution (cr, u) of (2.3)-(2.4) is the saddle point 
of the following energy 

E(t,v) = \ \\r\\ 2 + (divr,v) - (f,v), r e H(div;Q), v G L 2 (Q). 

Namely 

E(<t,u) — inf sup E(t,v). 

Similar result holds for the discrete solutions (cth, u h)- The lack of orthogonality is 
the main difficulty which complicates the convergence analysis of mixed finite element 
methods. 

We shall use the fact the error cr — cr H is orthogonal to the divergence free subspace of 
Tin to prove a quasi-orthogonality result. In the sequel we shall consider two conforming 
triangulations Th and Th which are nested in the sense that Th is a refinement of Th- 
Therefore the finite element space are nested i.e. Uh) C (E^, [/&.). 
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Lemma 3.1. Given an f 6 L 2 (Q) and two nested triangulation Th and Th, let 

(cr, u) = C' 1 /, (a h , u h ) = C^fh, {&h, u h ) = C^fn, and (a H , u H ) = C H l f H - 
Then 

(cr - cr h , cr h - cr H ) = 0. (3.1) 
Proof. Since cr h — cth & ^h, by (2.5)-(2.6), we have 

(cr - cr h , cr h - cr H ) = (u - u h , div (cr h - (T H )) = (u - u h , f H - Jh) = 0. 

□ 

To prove quasi-orthogonality, we need the following discrete stability result 

||cr h - a h \\ < vG)°sc(A,7h-), (3.2) 

where the constant C depends only on the shape regularity of Th- We shall leave the 
proof of (3.2) to the next section and use it to derive the quasi-orthogonality result. 

Theorem 3.2. Given an f 6 L 2 (fl) and two nested triangulations Th and Th, let cr = 
Z^f, °h = £ h X fh, andcr H = C^fn- Then 

{cr - cr h , cr h - cr H ) < \/Cq \W - cr h \\osc{f h , T H ), (3.3) 
Thus for any 5 > 0, 

(1 - 5)\\cr - cr h \\ 2 < \\cr - cr H \\ 2 - \\<T h - cr H \\ 2 + ^osc 2 (A, Th), (3.4) 

o 

and in particular when osc(//j, Th) = 0, 

II Il2 || ||2 || ||2 c\ 

WC — <Th\\ = \\°' — (t h\\ — \\<Th — <Th\\ ■ W-~>) 

Proof. Let us introduce an intermediate solution cr h = C h x f H - By Lemma 3.1, (cr — 
cr h , cr h - cr H ) = 0. Thus 

(cr — cr h , cr h — cr H ) = (cr — cr h , cr h — cr h ) < \\cr — cr h \\\\cr h — cr h \\. 

(3.3) then follows from the inequality (3.2). 

By the trivial identity cr — cr H = cr — cr h + cr h — cr H , we have 

II c — c h\\ 2 = || cr — c/i|| 2 + \\cr h — cr H \\ 2 + 2(cr — cr h ,cr h — cr H ) 

When osc(/fe,Tff) = 0, by (3.3), (cr — cr h ,cr h — cr H ) = and thus (3.5) follows. In 
general, we use 

II cr — v h\\ 2 = \\& — Ch || 2 + || Ch — cr H \\ 2 + 2(cr — cr h , cr h — cr H ) 

> \\cr - cr h \\ 2 + \\cr h - cr H \\ 2 - 2 V / C^||cr - cr h || osc(/, Th) 

> \\cr h - <t h \\ 2 + (1 - 5)\\cr - cr h \\ 2 -osc 2 (f h ,T H ), 



to prove (3.4). In the last step, we have used the inequality 

2ab < 5a 2 + \b 2 , for any 5 > 0. 
o 

□ 

A similar quasi-orthogonality result was obtained by Carstensen and Hoppe [26] for 
the lowest order Raviart-Thomas spaces using a special relation to the non-conforming 
finite element. Such relation for high order elements and Brezzi-Douglas-Marini spaces 
are not easy to establish; see [3] and [31, 32, 33] for discussion on this relation. In 
contrast the approach we used here is more straight-forward. 
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Remark 3.3. The oscillation term osc(// l , Th) in (3.3) and (3.4) depends on both Th and 
Th- It can be changed to the quantity osc(/, Th) which only depends on Th- Indeed for 
each T G Th, we have 

\\fh - f H \\ ,T = \\Qh(I - Q H )S\Wt < 11/ - f H \\o,T, 

and thus osc(fh,Tn) < osc(f,Tn)- This change is important for the construction of 
convergent AMFEM by showing the reduction of osc(/, Th)- 

4. Discrete stability for perturbation of data 

In this section, we shall prove the discrete stability result. We begin with a stabil- 
ity result in the continuous case. Let u E Hq(£1) be the solution of the primal weak 
formulation (2.1) of Poisson equation. Then (— Vu, u) is the solution to the dual weak 
formulation (2.3)-(2.4). The stability result ||V«|| < ||/||-i is well-known in the litera- 
ture. The norm ||/||_i, however, is not easy to compute. Instead we shall make use of 
the oscillation of data to bound it. 

Theorem 4.1. Given a shape regular triangulation Th ofQ and f E L 2 (Q), let (cr, u) = 
Cr x f and (<r, u) = Cr^jn, respectively. Then there exists a constant Co depending only 
on the shape regularity o/Th such that 

||o--<t|| < VC~oosc(f,T H ). (4.1) 

Proof. By (2.3) and (2.4), we have 

||<r — <t|| 2 = (cr — or, cr — a) = (div(cr — cr), u — u) = (/ — fn, u — u). 

Let v be the solution of primal weak formulation of Poisson equation with data / — fn- 
Then v = u — u and — Vf = cr — a. Recall that Qh '■ L 2 (Vl) — > Uh is the L 2 projection 
into discontinuous polynomial spaces. So for each triangle T G Th, (f — fn, vh)t = 
for any vh G V p (T). Therefore 

11°" - = (/ - fn,v) 

= ^2(f ~ fu,v-Q H v) T 

<Vc~oY1 II^(/-/»)||o,t||Vv|| ,t 

< V^(E W-/H)fc) 1/2 ||o--Cr||. 

TeT H 

In the second step, we have used the error estimate 

\\v - Qhv\\ ,t < ^/Q)H T \\ Vu|| ,t, 

which can be easily proved by Bramble-Hilbert lemma and the scaling argument. The 
constant Co only depends on the shape regularity of Th- The desired result then follows 
by canceling one || cr — cr || . □ 

In the proof of Theorem 4.1, we use the local error estimate 

||it - Qhu\\o,t < VCo-ffrllVwIlo/r = \ZCoH T \\cr\\ T , 

for u G Hq(VI) and cr = — Vw. The main difficulty in the discrete case is that u h E U h 
Hq(Q). However we still have a similar localized error estimate for Uh — QnUh- 



CONVERGENCE AND OPTIMALITY OF ADAPTIVE MIXED FE METHODS 9 

Lemma 4.2. Let Th and Th be two nested triangulations, and let (crh,Uh) = C^ 1 fh. 
Then for any T G Th, we have 

IK - QnUh\\o,T < \/~CQH T \\(T h \\^ T . (4.2) 

The proof of this lemma is technical and postponed to the end of this section. We use 
it to prove the following theorem. 

Theorem 4.3. Let Th and Th be two nested conforming triangulations. Let &h = 1 fn 
and cr h = C h l fh. Then there exists a constant C , depending only on the shape regularity 
ofTu such that 

II &h - &h || < V~Co osc(f h ,T H ). (4.3) 

Proof. Recall that <Jh — &h satisfies the equation 

{(T h - &h,T h ) = (u h - u h ,divT h ), Vr h G T, h (4.4) 

(div (<r h - a- h ),v h ) = (f h - fH,Vh), Wv h eU h . (4.5) 

We then choose r h = cr h — d"h in (4.4) and v h = Uh — Uh in (4.5) to obtain 

\\cr h - cr h \\ 2 = (u h - u h ,div((Th - &h)) = (v h , f h - fn) = (v h - Qnv h , fh - /#)• 

In the third step, we use the fact fjj = Qnf = Qnfh since Th and Th are nested. Thanks 
to (4.2), we have 

\Wh - °h\\ 2 = ^2 ( Vh ~ Qwh, fh - Jh)t 

TeT H 

< V~Q) H T \\f h — /h^tHc h — 07i||o,t 

Ter H 

< yCq I ^ H^Wfh - JhWt I \Wh - &h\\- 

\Te.T H ) 

Canceling one \\cr h — cr^||, we get the desired result. □ 

In the rest of this section, we shall prove Lemma 4.2. It is a modification of arguments 
in [4] from quasi-uniform grids to adaptive grids. The first ingredient is the existence of 
a continuous right inverse of the divergence as an operator from H" (fi) into the space 
Lg(fl) := {v G L\Q) : /„« = 0.} 

Lemma 4.4. Given a function f G Lj^fi), there exists a function r G H q(Q) such that 

divr = / and ||r||i < C\ 



The proof of this lemma for smooth or convex domains Vl is pretty easy. One can solve 
the Poisson equation with Neumann boundary condition 

A0 = / in Q, ^ = on dSl. 
on 

The condition / G L%(n) ensures the existence of the solution. Then we let r = grad0 
and modify the tangent component of r to be zero [22]. See also [7, 37] for a detailed 
proof on non-convex and general Lipschitz domains. 

The second ingredient is an interpolation operator Uh ■ H (fl) — > with the fol- 
lowing nice properties. 

Lemma 4.5. There exists an interpolation operator 11^ : H (T) — > such that 
(1) Q h div r = divn h r, Vr G H\n); 
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(2) there exists a constant C depending only on the shape regularity ofTh such that 

\\r - U h r\\ T < Ch T \\r\\ hT , VT G T h , Vr G H\n); 

(3) for any T G Th ifT G Hq(T), then HhT~\dT = 0. 

For the detailed construction of such interpolation operator and proof of these proper- 
ties, we refer to [41] and [6]. 

Proof of Lemma 4.2 We first note that Uh — QnUh = (Qh — Qh)uh since QhUh = Uh- For 
any T G Th, by the definition of L 2 projection Q H , we have, J T (Qh — Qii)uh = i.e. 
(Qh — Qh)uh G Lq(T). We thus can apply Lemma 4.4 to find a function r G H^(T) 
such that 

div T = (Q h - Q H )u h , in T and \\r\\ hT < C\\(Q h - Q H )u h \\ 0iT . 
We extend r to H 1 ^) by zero. Note that 

(n A - n H )r G S ft , and supp(n fe - U H )r C T. (4.6) 
With such r, we have 

IKQ/i - <2ff)wh||o,T = ((^ ~ Qn)u h , div r) T = Oh, (Qh - <5h) divr) T - 
Then using the commuting property (Lemma 4.5 (1)) and the locality of r, we have 

Oh, (Qh - Qh) div t) t = Oh, (Qh ~ Qh) div r) n = Oh, div(n ?i - )r) n . 

Now we shall use the fact (a h , u h ) is the solution of (2.3) and (2.4) and, again, the locality 
of r to get 

Oh, div(IT^ - U h )t)q = (cr h , (U h - U h )t) q = (cr h , (U h - U h )t) t . 
Using the approximation property of (Lemma 4.5 (2)), we get 

(cr h , (LLh - II h )t)t < ||<t/ 1 ||o,t(||t - U h r\\ ,T + \\r - n^r|| ,r) 

< Cf/ r ||<T/ 1 || ,T||''"||l ) T- 

So we have 

\\(Qh - QH)u h \\l T < Cflrlkftllo.rllTlli.T < CH T \\(T h \\ T \\(Q h - Q H )u h \\o }T - 
Canceling one || (Qh — Qn)uh\\T, we obtain the desired result. □ 

5. A Posterior Error Estimate for Mixed Finite Element Methods 

In this section we shall follow Alonso [1] to present a posteriori error estimate for 
mixed finite element methods. Other a posteriori error estimators for the mixed finite 
element methods can be found at [25, 62, 43, 39, 44, 45]. Our analysis could be adapted 
to these error estimators also. 

5.1. A posteriori error estimator and existing results. Let us begin with the definition 
of the error estimator. For any edge E G £h, we shall fix an unit tangent vector for 
E. We denote the patch of E consisting of triangles sharing E by Vt E . 

Definition 5.1. Given a triangulation l~H,for an E G £h and E ^ dfl, let Qe — T U T. 
For any cth G E#, we define the jump ofcrn across edge E as 

Je(cth) = \c h ■ ts\ '■= <J h\t ■ tE — cr h\t ' (5-1) 

If E G £h n dVL, we define Je(&h) = &h • tE- The edge error estimator is defined as 

V 2 e (ct h ) = \\Hrot(T H \\l nE + \\H^ 2 J e ((t h )\\Ie. 
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For a subset Th ^ £h> we define 

V 2 (<t h ,J 7 h) •■= ^e^h)- 

The error estimator tie((Th) is continuous with respect to <jh in L 2 -norm. Namely we 
have the following inequality. 

Lemma 5.2. Given an f G L 2 (Q) and a shape regular triangulation Th, let <r#, t# G 
There exists constant (3 such that 

(3 \r] 2 (cr H ,£ H ) - t] 2 (t h ,£h)\ < \\cr H - r H \\ 2 - (5.2) 

Proof. It can be easily proved by the triangle inequality and inverse inequality. □ 

We shall recall Alonso's results below and prove a discrete upper bound later. Since 
the data / is not included in the definition of our error estimator t)e, the upper bound 
contains an additional data oscillation term which is different from the standard one in 
[61]. 

Theorem 5.3 (Upper bound). Given an f G L 2 (Q) and a shape regular triangulation 
Th, let <t = C~ l f and cr H = L H X f H . There exist constants C and G\ depending only 
on the shape regularity ofTn such that 

Ik - (T H f < C iV 2 (<t Hi £ h ) + C osc 2 (f,T H ). (5.3) 

Theorem 5.4 (Lower bound). Given an f G L 2 (Q) and a shape regular triangulation 
Th, let (T = C~ l f and <th = C H l fn- There exists constant C% depending only on the 
shape regularity of Th such that 

C 2 v 2 (<r H ,E H ) < \\a-(T H \\ 2 , (5.4) 

for Raviart-Thomas spaces. 

For Brezzi-Douglas-Marini spaces, (5.4) holds when osc(/, Th) = 0. 

When osc(/, Th) = 0, (5.3) and (5.4) implies that C 2 /Ci < 1. This ratio is a measure 
of the precision of the indicator. 

5.2. Discrete upper bound. We shall give a discrete version of the upper bound (5.3). 
The main tool is the discrete Helmholtz decomposition. 
Given a shape regular triangulation Th, let 

SI = {4> h G C(U) : ij; h \ T G V P (T), VT G %} 

denote the standard continuous and piecewise polynomial finite element spaces of H 1 (fl) . 
To introduce the discrete Helmholtz decomposition, we define the dual operator operator 
of div : T, h U h . 

Definition 5.5. We define gr&d h : Uh i-> (S^)* by 

(gr&d h v h , r h ) = -{v h ,divT h ), Vr h G S fc . 

We emphasis that grad h is not simply the restriction of grad to Uh since Uh is not a 
subspace of H 1 ^). The following discrete Helmholtz decomposition is well known in 
the literature; See, for example, [38, 19, 5, 15]. 

Theorem 5.6 (Discrete Helmholtz Decomposition in IR 2 ). Given a triangulation Th, for 
p-th order Raviart-Thomas finite element spaces (S^, U%), we have the following orthog- 
onal (with respect to L 2 inner product) decomposition 

S^ = curl(^ +1 )©grad h (^). 
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For Brezzi-Douglas-Marini finite element spaces {Ti p h ,Uf l ), we have the following or- 
thogonal (with respect to L 2 inner product) decomposition 

E£ = curlOST 1 ) ©grader 1 )- 
We are in the position to present a discrete version of the upper bound. 

Theorem 5.7. Let Th and Th be two nested conforming triangulations. Let cr h = fh 
and (Th = fu, and let Th = {E G £ H : E ^ £ h }. Then there exist constants 
depending only on the shape regularity of Th such that 

\Wh - °h\\ 2 < C iV 2 (cth,Th) + C osc 2 (f h ,T H ) (5.5) 

and 

4F H <mT h -4T H ). (5.6) 
Proof. The inequality (5.6) follows from 

#T H < #£ h - #£ H < 3(#T, - #T H ). 

To prove (5.5), again we introduce the intermediate solution a h = Cf^fn- By the 
discrete Helmholtz decomposition, we have 

cr h -°H = grad h <\> h + curl^, 

where <p h G U^,if>h € S'f^ 1 for Raviart-Thomas spaces, and 4>h G Uf~ l ~,if>h € , for 
Brezzi-Douglas-Marini spaces. The decomposition is L 2 -orthogonal i.e. 

\\&h - cr H || 2 = || grad^f + || curl^f. (5.7) 

In two dimensions, || curled = || grad^Vill and thus (5.7) implies that 

\if>h\i < \\o-h ~ &h\\- (5.8) 

Since 

{& h - cr H ,grad h v h ) = (div (cr h - cr H ),v h ) = (f H - f H ,v h - v H ) = 0, 
we have 

\\<Th ~ vhW 2 = {o-h - cr H ,gTad h 4> h ) + (cr h - <t h , curl if) h ) = {& h - a H , curl if> h )- 

Since div curl if) = 0, (cr ft , curl^/J = and (cr H , curl = for any ip H G . 

Choosing ip H = ^H^h using some local quasi-interpolation, for example the Scott- 
Zhang quasi-interpolation [57], T H '■ S'h +1 h> S 1 ^ 1 such that 

1 /2 

\\lf) h -TBlf>h\\0,E < CHjj, \lf)h\l,Sl E and Wh -ZHlf>h\\Q,T < C H T \if>h\i,n T , 

where Qt = {Th C Th '■ Th flT ^ 0}. Furthermore the quasi-interpolation Th is 
local in the sense that if T G Th H Th or E G £ h D £h (i-e- T or E is not refined), then 
(iph — 1-Hif>h)\T = or (iph — ThiPh)\e = 0, respectively. With such choice of Th and 
ipH, we have 

W&h ~ Off || 2 = {<Th ~ <T H ,cur\ip h ) = (-cr H , curl (^ft - iPh)) 

TeT H E£dT 

< ^2l(T H ]\\^h-lpH\\o,E+ ^2 ll rotcr ffllll^ - 1pH\\o,T, 
E&£ h TeT H 

< Cr](a- H ,T H )\if>h\i < C^ifTmF^W&h - <th\\- 
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Canceling one \\cr h — cr H \\, we get 

\\&h-<rH\\<Cri(tr H ,F H ). (5.9) 
Now we write cr h — <r H = <r h — cr^ + cr fe — a h and note that 

{(T h - a h , cr h - (T H ) = (u h - u h , div(dr h - cr H ) = (u h - u h , f H - f H ) = 0. 
Combining (5.9) and (3.2), we then have 

\Wh - (?h\\ 2 = W&h - cr H \\ 2 + \\a h - cr h \\ 2 < Cirj 2 (cr H ,J z H ) + C osc 2 (f h , Th)- 

□ 

6. Convergence and Optimality of AMFEM 

In this section we shall present our algorithms and prove their convergence and opti- 
mality. It is adapted from the literature [48, 59, 36, 49, 50]. For the completeness we 
include them here and prove some important technical results. 

We first present our algorithms. It mainly follows from the algorithm proposed in [50] . 
The difference is that we do not impose an interior point property in the refinement step. 

Let % be a initial shape regular triangulation, a right hand side / G L 2 (Vt), a tolerance 
e, and < 9,9, fi < 1 three parameters. Thereafter we replace the subscript h by an 
iteration counter called k. For a marked edge set Aik, we denote by Vt Mk = VJe^mS^e- 

[T n , cr JV ]=AMFEM(r , /, e, 9, 9, //) 

rj = e, k = 

WHILE 7] > e, DO 

Solve (2.5)-(2.6) on % to get the solution cr k . 
Compute the error estimator rj = r/(cr fc , £ k ). 
Mark the minimal edge set M. k such that 

r, 2 {(T k ,M k ) > 9r ] 2 ((T k ,S k ). (6.1) 

If osc(/, Tk) > osc(/, 7o)/i fc , enlarge M k such that 

osc(/,OmJ > 9osc(f,T k ). (6.2) 

Refine each triangle r G &M k by the newest vertex bisection to get T+i- 
fc = Jfe + l. 
END WHILE 

Tn = T- 
END AMFEM 

6.1. Convergence of AMFEM. We shall prove the algorithm AMFEM will terminate 
in finite steps by showing the reduction of the sum of the error and the error estimator. 

We first summarize the main ingredients in the following lemma with the following 
short notation: 

e k = || cr - (T k \\ 2 ,E k = \\cr k +i - cr k \\ 2 ,o k = osc 2 (f,T k ), and i] k = r] 2 ((T k ,£ k ). 



Lemma 6.1. One has the following inequalities 

C 

(1 - 5) e k+1 < e k - E k + — o k , for any 5 > (6.3) 

o 

/3r) k+1 </3{l-^0)r] k + E k , (6.4) 

e k < G x r] k + Co o k - (6.5) 
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Proof. (6.3) is the quasi-orthogonality (3.4) established in Theorem 3.2 and Remark 3.3. 
(6.5) is the upper bound (5.3) in Theorem 5.7. We only need to prove (6.4). By the 
continuity of the error estimator (5.2), we have 

f3r] 2 ((T k+1 ,£ k+1 ) < /V(°"fc> £k+i) + E k . (6.6) 

Let jVfe+i = £k+i\£k be the new edges in T k +i and M k C £ k be the refined edge in %. It 
is obvious that £ k \M k = £ k+ i\M k+ i. For an edge E G A4+i, if it is an interior edge of 
some triangle T E %, then Js(cr fe ) = since cr k is a polynomial in T. For other edges, 
it is at least half of some edge in Ai k and thus 

7] 2 {(T k ,Af k+1 ) < -r] 2 {(T k ,M k ). (6.7) 

Since some edges are refined for the conformity of triangulation, Ai k C M. k . By the 
marking strategy (6.1), we have 

rf(<r kl M k ) > V 2 {(T k ,M k ) > 6 V 2 {(T k ,£ k ). (6.8) 

Combining (6.7) and (6.8), we get 

7] 2 ((T k ,£ k+1 ) = i] 2 (cr k ,J\f k+1 ) +r] 2 (cr k ,£ k+1 \J\f k+1 ) 

< ^n 2 (<T k ,M k ) +7] 2 {(T k ,£ k \M k ) 

< -^(vk, M k ) + T] 2 ((T k , £ k ) 

< (l-i0)77vk,£fc)- 

Substituting to (6.6) we then get (6.4). □ 
Theorem 6.2. When 

0<5<min{^-M}, (6.9) 
there exists a G (0,1) and C$ such that 

(1 - 8)e k+ i + pr] h+1 < a[(l - 5)e k + pr} k ] +C s o k . (6.10) 
Proof. First (6.3) + (6.4) gives 

1 C 

(1 - 5)e k+1 + (3r] k+1 < e k + (3(1 - -6)rj k + -z- o fc . 

2 o 

Then we separate e k and use (6.4) to bound 

e k = a(l - 5)e k + [1 - a(l - 5)]e fc 

< a(l - 5)e fc + [1 - a(l - <S)](Ci7fo + C o fc ). 

Therefore we obtain 

(1 - <5)e fc+1 + < a |(1 - <5)e fe + ^ ~ + Q o*. 



Now we choose a such that 

[1 - a(l - 8)}C X 



a 

i.e. 



Ci + (1 - \9)P Cx + p-h 

a 



Cx(l-5)+p Cx + p-CxS' 
By the requirement of 5 (6.9), we conclude a G (0, 1). □ 
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Theorem 6.3. Let cr k be the solution obtained in the k-th loop in the algorithm AMFEM, 
then for any < 5 < minl^-^, 1}, there exist positive constants C$ and < 75 < 1 
depending only on given data and the initial grid such that, 

(1 - 5)11(7 - <r k \\ 2 + 0r) 2 (<r k , %) < C 5l k Sl 

and thus the algorithm AMFEM will terminate infinite steps. 

Proof. The proof is identical to that of Theorem 4.7 in [50] using (6.10). □ 

6.2. Optimality of AMFEM. Let 7o be an initial quasi-uniform triangulation with #7o > 
2 and V N be the set of all triangulations T which is refined from % and #7" < N. For a 
given triangulation T, the solution of the mixed finite element approximation of Poisson 
equation will be denoted by cr T . We define 

A s = {cr G £ : \\cr\\ A s < 00, with \\cr\\ A s = sup (N s inf \\cr - <r r ||)}. 

An adaptive finite element method realizes optimal convergence rates if whenever cr e 
A s , it produces approximation cr N with respect to triangulations Tn elements such that 
\\cr-cr N \\ <C{#T N )- S . 

For simplicity, we consider the following algorithm which separates the reduction of 
data oscillation and error. 

(1) [T H , fu] = APPROX (/, To, e/2) 

(2) [cr N , T N ] = AMFEM (Tff, f H , e/2, 6, 0, 1) 

The advantage of separating data error and discretization error is that in the second 
step, data oscillation is always zero since the input data /# is piecewise polynomial in 
the initial grid Th for AMFEM. In this case, we also list all ingredients needed for the 
optimality of adaptive procedure. 

(1) Orthogonality: 

II I|2 || ||2 || |2 

W — Cfc+i|| = W — <f fc|| — Wk+l — 0"fc|| 

(2) Discrete upper bound: 

||cr fc +i - cr k \\ 2 < Cirfi^cTk^k) and #J" fc < 3(#7^ + i - #7^). 

(3) Lower bound: 

C 2 rj 2 {(T k ,£ k ) < \\cr - cr k \\ 2 . 

Theorem 6.4. Let [cr N , T N ] = AMFEM(Th, f H , e, 6, 0, (i), and cr = L~ x f H - If a G A s 
and < 6 < C 2 /Ci, then for any e > 0, AMFEM will terminated infinite steps and 

\\cr - cr N \\ <e, and #T N - #T < C\\<r\H' e~ y ' ■ (6.11) 

Proof. It is identical to the proof of Theorem 5.3 in [59] using three ingredients listed 
above. □ 

Theorem 6.5. For any f e L 2 (Q), a shape regular triangulation % and e > 0. Let cr = 
C- l fand [cr N , T N ] = AMFEM (T H , f H , e/2, 0, 1) where [Th, f H ] = APPROX (/, %, e/2). 
IfcrEA s and f G A s a , then 

Ik - cr N \\ < C(\\cr\\ A s + ||/|U s )(#71v " #T )- S . 
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Proof. Let cr = C }h- By Theorem 4.1 and 2.2, we have 

\\a-a\\<e/2, and #T H - #T < CWffp- 1 ^. 
It is easy to show, by the definition of A s , if cr e .4 s , then cr e .4 s and 

II^IU" < IMU* + ll/IUi 

We then apply Theorem 6.4 to or to obtain 



(6.12) 



o- - (t n \\ < e/2 and #T W - < C||a-||^- 1/s 



(6.13) 




e < C(#71v - #7 )- s (|| C r|U s + H/lUg). 



The desired result then follows. 



□ 



7. Conclusion and future work 



In this paper, we have designed and analyzed convergent adaptive mixed finite ele- 
ment methods with optimal complexity for arbitrary order Raviart-Thomas and Brezzi- 
Douglas-Marini elements. Although the results are presented in two dimensions, most of 
them are dimensional independent. For example, the discrete stability result, Theorem 
4. 1 , holds in arbitrary dimensions without any modification of the proof. 

The proof for the upper bound of the error estimator (Theorem 5.3 and 5.7), however, 
cannot be generalized to three dimensions in a straightforward way. In the proof, we use 
a special fact that in two dimensions, if (curl) is as smooth as H 1 since in two dimensions 
curl operator is just a rotation of gradient operator. To overcome this difficulty, we need 
to use a regular decomposition instead of Helmholtz decomposition. Note that discrete 
regular decomposition for corresponding finite element spaces is developed recently by 
Hiptmair and Xu [42]. We could use these techniques to prove the convergence and 
optimality of adaptive mixed finite element methods in three and higher dimensions. 

Acknowledgement. The authors would like to thank Dr. Guzman for the discussion on 
the simplification of the proof of the discrete stability result and Prof. Nochetto for the 
simplification of convergence analysis without using discrete lower bound. 
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